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ABSTRACT 

We present the Superparticle Model/ Algorithm for Collisions in Kuiper belts and debris disks 
(SMACK), a new method for simultaneously modeling, in 3-D, the collisional and dynamical evo- 
lution of planetesimals in a debris disk with planets. SMACK can simulate azimuthal asymmetries 
and how these asymmetries evolve over time. We show that SMACK is stable to numerical viscosity 
and numerical heating over 10' yr, and that it can reproduce analytic models of disk evolution. We 
use SMACK to model the evolution of a debris ring containing a planet on an eccentric orbit. Dif- 
ferential precession creates a spiral structure as the ring evolves, but collisions subsequently break up 
the spiral, leaving a narrower eccentric ring. 

Subject headings: Celestial mechanics — circumstellar matter interplanetary medium - methods: 
numerical — planetdisk interactions — planetary systems 


1. INTRODUCTION 

Spatially resolved debris disk images from optical and 
infrared observatories show spectacular patterns and 
sub-structures including eccentric rings (e.g. Kalas et al. 
2005; Schneider et al. 2009; Krist et al. 2012; Boley et al. 
2012), warps or sub-disks (e.g. Golimowski et al. 2006; 
Krist et al. 2005), and various other morphologies (e.g. 
Hines et al. 2007; Kalas et al. 2007). Undetected exoplan- 
ets could create many of these features via gravitational 
perturbations. Many authors have analyzed resolved im- 
ages of debris disks to predict the presence of exoplanets 
and constrain their locations, orbits, and physical prop- 
erties (e.g. Wyatt et al. 1999; Greaves et al. 2005; Quillen 
2006; Stark & Kuchner 2008). In the last few years, di- 
rect images of exoplanets combined with numerical mod- 
els (e.g. Chiang et al. 2009; Lagrange et al. 2010) have 
demonstrated the power and necessity of this approach. 
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However, our ability to use debris disk asymmetries as 
signposts of planets is limited by our ability to model 
debris collisions. The collisional lifetime of a dust grain 
with orbital period t per in a disk of optical depth r e ff 
is given by t co i = t per /4:TTT e f / (Wyatt 2008). In many 
debris disk systems, this collisional timescale is shorter 
than the timescales for dynamical sculpting by planets, 
so collisions can have a significant effect on the dynamics 
of the disk compared to the gravity of the planet. 

For example, Chiang et al. (2009) estimated the mass 
of a planet sculpting the Fomalhaut ring to be M p i ss 
0.5 Mj and the ring to have an optical depth of r e // ~ 
10 -3 . The secular dynamical perturbations from a planet 
of mass M p i act on a timescale of t sec ss t per M*/M p i 
(Murray & Dermott 1999) but would have a collisional 
lifetime of t co i ~ 3.8 x 10 4 yr, so a planetesimal orbiting 
Fomalhaut at 140 AU would experience secular effects on 
a timescale of t sec w4x 10 6 yr. Resonant effects, which 
could produce azimuthally asymmetric structures, can 
act on comparable timescales as well (Kuchner & Holman 
2003). In the Fomalhaut disk, the resonant timescale is 
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Ues « t per (M r / M^) 1 / 2 W 7 X 10 4 yr. Models of secular 
and resonant phenomena in such systems must therefore 
take collisions into account. 

There have been many attempts to wed collisional and 
dynamical effects in a debris disk model. Early models 
adapted N-body simulations with routines that spawned 
new code bodies at every collision. For example, Beauge 
& Aarseth (1990) represented the products of each frag- 
mentation with five daughter bodies. Grigorieva et al. 
(2007) represented fragmentation products by introduc- 
ing one daughter body per decade in mass. However, 
in this kind of algorithm, the number of code particles 
quickly increased to unmanageable numbers, so these 
models could not be run for very many orbits. 

Other modelers took the opposite approach: they be- 
gan with particle-in-a-box codes tracking averaged dy- 
namical quantities and added spatial resolution and other 
refinements to increase the dynamical fidelity. For ex- 
ample, Kenyon & Bromley (2006) utilized a disk di- 
vided into a series of rings, each of which contained a 
particle-in-a-box calculation. The ACE code (Krivov 
et al. 2005; Lohne et al. 2008; Vitense et al. 2012) utilizes 
the Boltzmann-Schmoluchowski equation to model colli- 
sional and dynamical evolution in debris disks with dy- 
namics averaged over the angular orbital elements. These 
kinds of codes can simulate 10 9 yr of disk evolution on 
existing computers (e.g. Vitense et al. 2012), but cannot 
be used to model disk asymmetries in three dimensions. 

Still other methods employ fully three-dimensional dy- 
namics, but only aim to model systems that are in steady 
state, i.e., the sources and sinks of dust grains roughly 
balance one another, as they might in a system that has 
already evolved for many collision times. One method 
that makes this approximation is the collisional groom- 
ing algorithm (e.g. Stark & Kuchner 2009; Kuchner & 
Stark 2010). Another is DyCoSS (Thebault 2012). 

A new model, LIPAD (Levison et al. 2012), uses a su- 
perparticle method coupled with full N-body integrators 
to simulate collisional and dynamical evolution, while 
keeping the number of superparticles roughly constant. 
Using this approach avoids the steady-state approxima- 
tion. We will discuss superparticle methods in Section 

2 . 1 . 

But despite all this work, no published colli- 
sional/dynamic disk models quite met our desires for in- 
terpreting images of debris disks. We wanted a model 
that could 

• Track collisional and dynamical evolution of debris 
disks in 3-D to model asymmetries created by plan- 
ets like warps and eccentric rings, and 

• Run stably for 10 7 or more years of disk evolution 
in a feasible number of CPU cycles. 

The above algorithms did not meet these requirements. 
For example, the LIPAD code uses a scale-height approx- 
imation for the disk’s vertical structure, and the longest 
published LIPAD simulations ran for only 10 4 yr. 

Therefore, we have developed a new tool for 3-D model- 
ing of collisional planetesimal populations in debris disks. 
Our tool, the Superparticle-Method Algorithm for Colli- 
sions in Kuiper belts and debris disks (SMACK), uses a 
superparticle approximation to simultaneously track the 


N-body dynamics and collisional evolution of the bod- 
ies that produce the dust we observe and calculate the 
dust production rates. We have designed SMACK with 
the ultimate goal of deriving improved estimates of the 
masses and orbital parameters of exoplanets in debris 
disks using high spatial resolution images, e.g., from the 
Atacama Large Millimeter Array (ALMA). This paper 
describes the SMACK algorithm, and presents a basic 
SMACK model for an eccentric debris ring. 

2. THE MODELING TOOL 

2.1. Superparticles 

Our full numerical model uses an N-body integrator, 
REBOUND, to solve the equations of motion of the plan- 
etesimals and detect collisions, combined with a collision 
resolution algorithm, SMACK, to calculate the effects 
of collisions on the velocities and size distributions of 
the planetesimals. REBOUND (Rein & Liu 2012) is a 
multi-purpose code originally designed for studying col- 
lisional dynamics in planetary rings, freely available un- 
der an open-source license from https://github.com/ 
hannorein/rebound. While the Grigorieva et al. (2007) 
model and LIPAD both detect collisions using a two- 
dimensional grid, REBOUND detects collisions in 3-D, 
without using a grid; it contains a Barnes-Hut tree mod- 
ule to calculate self-gravity and detect collisions, paral- 
lelized with MPI and OpenMP. 

In the standard version of REBOUND, collisions are 
resolved using an instantaneous collision model with a 
normal coefficient of restitution. This approximation 
handles a variety of problems in planetary ring dynam- 
ics. Debris disks contain a very different physical regime 
than planetary rings; collision speeds in debris disks are 
much higher than in rings (km s _1 vs. mm s _1 ) and 
optical depths are lower (~ 10 -4 vs ~ 1). So rather 
than allowing each REBOUND code body to represent 
one planetesimal, we use each body as either a planet 
with mass, or as a massless “superparticle” representing 
a group of planetesimals on similar trajectories. 

Superparticle methods have already been used exten- 
sively to model planetesimal formation within gas disks 
(e.g Michikoshi et al. 2009; Rein et al. 2010; Zsom & 
Dullemond 2008; Johansen et al. 2012; Charnoz & Tail- 
lifet 2012). They have also been applied to model 
planetesimals in the solar system and in debris disks. 
Charnoz & Morbidelli (2003) used a superparticle ap- 
proach to study how the dynamics of particles ejected 
from the Jupiter-Saturn region affected their size distri- 
butions, though they did not include the feedback from 
the collisions on the particle dynamics. Grigorieva et al. 
(2007) used cylindrical superparticles 5 AU in diameter 
fixed to the disk midplane to model collisional avalanches 
in debris disks over spans of ~ 40 orbital periods. The 
LIPAD model (Levison et al. 2012) uses a superparticle 
approach; these authors refer to their superparticles as 
“tracers” . 

2.2. SMACK 

In SMACK, the superparticles represent collections of 
planetesimals with a range of sizes, as in Charnoz & Mor- 
bidelli (2003). However, in Charnoz & Morbidelli (2003), 
the evolution of the size distribution of each superparticle 
is calculated only after the entire dynamical integration 


Collisions in Debris Disks 


3 


of the superparticles is complete. In contrast, SMACK 
calculates the size evolution at each timestep of the N- 
body integration, allowing us to model the feedback be- 
tween the collisions and the dynamics of the superpar- 
ticles. Each superparticle in SMACK is characterized 
by an incremental size distribution with logarithmic size 
bins. 

In general, two colliding, fragmenting bodies produce 
a spray of daughter particles with different trajectories. 
Simulating this distribution of trajectories in great de- 
tail would require SMACK to create new superparticles 
with each collision, quickly increasing the number of bod- 
ies tracked by REBOUND with every collision. Instead, 
SMACK approximates the outcomes of collisions while 
keeping the number of integrator bodies constant. When 
REBOUND detects an encounter between two superpar- 
ticles, it passes the velocities and size distributions of 
the overlapping superparticles to SMACK. SMACK re- 
turns two superparticles with different velocities and size 
distributions and REBOUND continues its dynamical in- 
tegrations using these modified superparticles. 

The essence of SMACK is simple. In SMACK, the 
fragments produced by collisions are swapped between 
superparticles, so that the new size distributions n a and 
n b are given by 

n a (i) = n A {i) - P A {i) + F B {i) (1) 

n b {i) = n B (i) - Pb (*) + F A (i), (2) 

where P A (i ) is the number of parent body particles in 
size bin i in superparticle A that are lost due to colli- 
sions, F b {%) is the number of daughter particles in size 
bin i produced by colliding particles in superparticle R, 
and n A and n B are the size distributions of the parent 
superparticles. The detailed forms of P(i) and F(i) used 
in SMACK are given in Section 2.3. 

This swapping of fragments is a first-order approxi- 
mation of the velocity distribution of fragments in a 
planetesimal collision. In a collision between two real 
planetesimals, the fragments are produced in roughly the 
center-of-mass frame. But an encounter between two su- 
perparticles is more complicated; the center of mass of 
the superparticles is not the same, in general, as the cen- 
ter of mass of any pair of planetesimals represented by 
the superparticles. If the mass ratio of the parent bodies 
is higher (lower) than that of the superparticles, the frag- 
ments should be launched in a direction skewed toward 
that of the more (less) massive parent bodies. The swap- 
ping described above crudely approximates this physics. 

After swapping the daughter planetesimals, SMACK 
calculates the kinetic energy lost in the collisions and 
corrects the superparticle velocities to reflect this loss 
and to conserve momentum. Let v A and v B be the mag- 
nitudes of the velocities of the parent superparticles in 
the center-of-momentum frame and m A and m B be the 
total masses of the parent superparticles. Some fraction 
of the kinetic energies of the parent superparticles will 
be lost to planetesimal collisions. This fraction depends 
on the fraction of planetesimals that experience collisions 
and the amount of kinetic energy lost by each planetes- 
imal in a collision. Using the energy loss fraction and 
the collision rates for each pair of size bins, SMACK cal- 
culates E a and E b , the kinetic energy lost to collisions 
in superparticles A and B, respectively (as described in 


Section 2.3). Then the new superparticle velocities must 
satisfy the energy conservation law, 

K = ^m a vl + ^m b v b = : m A v] \ + \m B v | - E A - E B , 

( 3 ) 

where K is the total kinetic energy of the two superparti- 
cles, and m a and m b are the new total masses of each su- 
perparticle, calculated from the new post-encounter size 
distributions given by equations (1) and (2). The veloc- 
ities must also satisfy the momentum conservation law, 

m a v a + m b v b = 0 (4) 

The velocities that solve Equations (3) and (4) are 


2 m b K 

m a (m a + m b ) () 

™a 

v b = v a . (6) 

m b 

Equations (1), (2), (5), and (6) give the new size distri- 
butions and velocities of the superparticles, and define 
the essential SMACK algorithm. 

When a superparticle encounter yields no planetesi- 
mal collisions (i.e., all the planetesimals pass through 
the encounter unaffected), the algorithm gives exact re- 
sults. When all the planetesimals in each superparticle 
collide, the algorithm produces an outcome that is a good 
approximation for the dominant size bins. When there 
is a mix of collisions and pass-throughs, the algorithm 
compromises, making errors in the distribution of out- 
put velocities, but not in the total energy or angular 
momentum. 

2.3. Fragmentation 

For now, SMACK only models one type of collision 
outcomes: catastrophic collisions, defined as collisions 
in which the largest fragment is no larger than half the 
size of the target. Cratering collisions do not have a sig- 
nificant effect on the steady-state size distribution of a 
collisional cascade (Dohnanyi 1969), so we neglect them 
for now. Since we are modeling disks in which the im- 
pact velocities are high, we also ignore bouncing, stick- 
ing, and gravity between planetesimals. Future versions 
of SMACK may incorporate these effects. 

Consider two superparticles, A and B, that are found 
to overlap during a given timestep. The number of plan- 
etesimals from size bin i in superparticle A that collide 
with planetesimals in size bin j in superparticle B is 

c A (i,j) =n A (i)T B (i,j), (7) 

where n A (i) is the number of planetesimals in size bin i 
in superparticle A and r B (i,j) is the optical depth along 
the path of a planetesimal in size bin i passing through 
superparticle B for collisions with planetesimals of size 
j. SMACK estimates this optical depth as 

r s(bi) ~ n B (j)crijl/V, (8) 

where er^- is the combined collisional cross-section of a 
planetesimal in size bin i and a planetesimal in size bin 
j, l is the path length traveled by a planetesimal in su- 
perparticle A, and V is the volume of superparticle B. 


4 


Nesvold et al. 


The path length l is the distance traveled by superparti- 
cle A through the disk relative to the local Keplerian flow 
since its last encounter. SMACK estimates this distance 
as l = VABtenc , where vab is the magnitude of the rel- 
ative velocity of superparticles A and B, and t enc is the 
time since superparticle A’s last encounter. The super- 
particle volume V is the parameter used by REBOUND 
to determine when a collision occurs. The superparticle 
volume must be chosen carefully to minimize both com- 
putation time (see Section 2.4) and numerical heating 
(see Section 3.2). We run numerical heating tests before 
every new simulation to select the optimal superparti- 
cle size. The collisional cross-section is purely geometric 
because gravitational effects between planetesimals are 
ignored. The cross-section is given by 

°H = l(D(i) + D(j)) 2 , (9) 

where D{i) and D(j) are the diameters of planetesimals 
in size bins i and j, respectively. 

We aspire to model a system in which each planetes- 
imal involved in a fragmenting collision loses some frac- 
tion fxE of its kinetic energy. This fraction fxE is a 
parameter of the code. The total desired energy loss for 
size bin i in superparticle A during its encounter with 
superparticle B is calculated by SMACK using 

e a(i) =T,fKE^m{i)v A (i,j) 2 c A (i,j), ( 10 ) 

3 

where VA(i,j) is the magnitude of the velocity of super- 
particle A in the center-of-momentum frame of a plan- 
etesimal in size bin i in A and a planetesimal in size bin 
j in B. 

Some of the colliding planetesimals counted in equa- 
tion (7) may shatter, creating a distribution of smaller 
fragments. SMACK determines which colliding planetes- 
imals will fragment by comparing the collision energy of 
each pair of colliding planetesimals to the minimum ki- 
netic energy needed for a catastrophic collision. In the 
center-of-mass frame, the kinetic energy of two colliding 
planetesimals with masses m{i) and m(j) is 


where G is the gravitational constant, m(i) is the mass 
of planetesimals in size bin i, D(i) is the diameter of 
planetesimals in size bin i, and S is the impact strength 
of the planetesimals. We use }ke = 0.1 (Fujiwara 1982) 
and a size-independent impact strength of S = 3 x 10 6 J ■ 
m~ 3 (Greenberg et al. 1977). 

If Equation (12) holds true, the loss of parent bodies in 
size bin i due to fragmentation in collisions with size bin j 
is equal to the number of collisions between planetesimals 
in i and j: 

Pa{i, j) = c A (i, j )■ (14) 

If the inequality in Equation (12) does not hold, the plan- 
etesimal in size bin i will not shatter and will not be 
counted as loss, so Pa (i, j) = 0. The planetesimal in 
size bin j may fragment, however, depending on whether 
E m in{j) satisfies Equation (12); SMACK performs these 
calculations by looping through one index at a time. The 
total loss in size bin i in superparticle A is the sum of 
the losses due to each size bin in superparticle B: 

p a{i) = Y^ PA (iJ)- ( 15 ) 

j 

Collisions produce fragments in power law size distri- 
butions, which we represent with incremental logarith- 
mic bins. The individual daughter particle distribution 
resulting from fragmentation in size bin i is 

F A{iJ) = PA(i,j)nAD{i) a . (16) 

The fragment distribution index, a, is a parameter of 
the algorithm. The constant ka is calculated for ev- 
ery collision for each superparticle such that the largest 
fragment in each distribution is one-half the mass of the 
parent planetesimal that produces the fragments. Again, 
we find the total fragment gain from size bin i in super- 
particle A by summing over the size bins in superparticle 
B that collide with i: 

F A (i) = Y, F A(iJ)- (17) 
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2.4. The Smallest Planetesimals 


1 m(i)m(j) 2 
= 2m(i) + m{j) Vrel ' 


( 11 ) 


where v re i is the magnitude of the relative velocity of the 
two planetesimals. In laboratory experiments, Hartmann 
(1980) found that approximately half of this collisional 
kinetic energy is partitioned into each colliding body, re- 
gardless of the mass ratio of the two, so a planetesimal 
in size bin i will shatter in a collision with a planetesimal 
in size bin j if 

2 E col(hj, Vrel) A E min{l')-i ( 12 ) 


where E min (i ) is the minimum energy needed to shatter 
a planetesimal in size bin i. 

The minimum shattering energy, E m i n (i), is a combi- 
nation of the gravitational binding energy of the plan- 
etesimal and its internal impact strength. We use the 
minimum energy criteria derived by Durda (1993), 


Emin(i') — 


fKE 


0.822 


Gra(i) 2 


X7T SD(tf 

O 


(13) 


If the total optical depth for a given size bin is ever 
greater than one, i.e. , 

Tb(*) = ^2 T B {hj) > U (18) 

3 

the loss in that size bin given by equation (15) could be 
greater than the number of planetesimals in that bin. 
We generally try to avoid this situation by choosing the 
number and volume of the superparticles so that the su- 
perparticle encounter time is less than the collision time 
for the small grains. However, it sometimes occurs any- 
way as the disk evolves and more small planetesimals 
are produced. This situation tends to arise for the small- 
est planetesimals first, since they collide most frequently, 
and become more common with increasing superparticle 
encounter times and with increasing optical depths. 

To address this issue, we adopted a variable-timestep 
method. If the maximum optical depth for all the size 
bins in a superparticle is ever found to be greater than 1, 
SMACK divides the path length traveled by the super- 
particle since its last encounter into segments such that 
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the maximum optical depth in any size bin is < 1. Af- 
ter each segment, the optical depths are recalculated and 
the energy and planetesimal losses are calculated as de- 
scribed above. The relative velocity of the superparticles 
is updated after each segment to accurately reflect the 
energy lost to collisions, but the superparticle trajecto- 
ries are not changed until the final segment. 

This “small particle loop” allows us to run SMACK un- 
interrupted in regions of unexpectedly high density with 
fewer superparticles. However, it also introduces noise 
into the velocity evolution of the superparticles, as the 
trajectories are not updated during the loop. We there- 
fore check each astrophysical simulation by running it 
with a few different superparticle sizes to ensure that we 
get the same result, and that the small particle loop is 
not introducing excess noise. 

While REBOUND can be adapted to model small dust 
grains by adding addition forces to the integrator such as 
radiation pressure and Poynting-Robinson drag, SMACK 
cannot simultaneously model dust grains and planetesi- 
mals within the same superparticles if they are subject 
to different forces. We therefore model only grains larger 
than 1 mm, which do not experience significant radiative 
forces during their collisional lifetimes, even in a very 
sparse disk like the zodiacal cloud. The current version 
of SMACK thus cannot directly simulate disk images at 
optical wavelengths, but is meant for simulating high- 
resolution millimeter and sub-millimeter images from in- 
struments like ALMA. 

2.5. Normalization 

To compare our models with images and photometry 
of known debris disks, we need to know the face-on opti- 
cal depth, Tdi S k, of the models, and how it relates to rgp, 
the total cross section of the planetesimals in the super- 
particle divided by the cross section of the superparticle, 
7T r^. Equation (8) above implies that the density of an 
individual superparticle models the local disk density. So 
T disk is related to the individual optical depth of each su- 
perparticle by a linear filling factor, fsp, where 

Tdisk = fsP T SP- (19) 

The filling factor, fsp, is simply the average number 
of superparticles that a perpendicular path through the 
disk would intersect. In other words, fgp is the inverse of 
the fraction of the perpendicular path through the disk 
that would be contained within one superparticle. This 
factor is approximately 

( 20 ) 

where h is the full height of the planetesimal distribution 
and 4r s /3 is the average length of a chord through a 
spherical superparticle with radius r s . 

We calculate fsp before the simulation begins using 
equation (20). We generate 10 6 superparticles with the 
same orbital element distributions that we will use as the 
initial conditions. Then we choose a fiducial perpendic- 
ular path through the disk, and create a histogram of 
the positions of the superparticles along that path. We 
estimate h along that path by normalizing the histogram 
such that the maximum value is 1 and summing together 
all the bins. Knowing fgp allows us to set the total cross 


section of the planetesimals in the superparticles to yield 
an initial face-on optical depth of our choosing. 

The planetesimal size distributions in the superparti- 
cles begin at 1 mm, but dust grains in disks are ground 
down to ~ 1 /im in size before they are removed from 
the system by radiation pressure. To ensure that we are 
including the contribution of these smaller dust grains to 
the cross-sectional area of the planetesimals in the disk, 
we fit a power law to the size distributions of each super- 
particle, then extrapolate the size distributions down to 
I /rm. We use the extrapolated size distributions when 
calculating the cross-sectional area of the superparticles 
during normalization. 

3. NUMERICAL TESTS 

We performed a series of numerical tests on SMACK to 
validate the algorithm and identify sources of numerical 
noise. We used the Wisdom-Holman integrator (Wis- 
dom & Holman 1991) included in REBOUND, which 
closely follows the implementation of the SWIFT code 
(Levison & Duncan 1994). For collision detection, we 
selected REBOUND’S tree algorithm, which implements 
a nearest-neighbor search to find overlapping superparti- 
cles at each timestep. We ran REBOUND and SMACK 
on the NASA Center for Climate Simulation’s (NCCS) 
Discover cluster, using a hybrid OpenMP/MPI paral- 
lelization on 48 cores. 

In each simulation, we assume the planetesimals are 
spherical with a density of 3 g/cm 3 . We also assume that 
the power law distribution of the fragments (see equa- 
tion 16) has index a = —2.8. We also assume the star 
has mass 1 Mq and radius 1 Rq. REBOUND outputs 
the orbital elements, Cartesian coordinates, and size dis- 
tributions of each superparticle every output timestep, 
where the size of the output timestep is greater than the 
integrator timestep and is set by the user. For each sim- 
ulation, we use an integrator timestep of 1 yr and an 
output timestep such that REBOUND outputs 1000 to- 
tal data points for each superparticle. For example, we 
use the output timestep of 10 4 yr for a 10 7 yr simulation. 
We use open boundary conditions for our systems; if at 
any timestep a superparticle’s orbital elements place it 
outside a cube with a user-defined width h, ox centered 
on the star, the superparticle is removed from the sim- 
ulation. The system boundaries form a cube because 
REBOUND was originally developed for studying shear- 
ing boxes. If a superparticle collides with the star or any 
planet in the system, the superparticle is also removed. 
The initial conditions for each simulation presented in 
this section are shown in Table (1). The angular orbital 
elements (longitude of the ascending node f2, argument 
of pericenter w, and true anomaly /), are all distributed 
uniformly between 0 and 2-7t for each simulation. 

3.1. Size Distribution Evolution 

Dohnanyi (1969) studied collisional cascades analyti- 
cally assuming an infinite range of particle masses with 
each body experiencing a constant impact velocity, as- 
suming a mass-independent material strength. He found 
that the incremental mass distribution of such a colli- 
sional cascade at steady state can be described by a 
power law with index q = —1.833 (with linear mass bins). 
The corresponding incremental size distribution with log- 
arithmic bins is a power law with index p = —2.5 (Durda 
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1993). 

Dohnanyi (1969) also found that the index of the equi- 
librium size distribution is independent of the fragment 
size distribution index a. Durda (1993) confirmed this in 
his collisional model, and noted that steeper values of a 
corresponded to a faster convergence of the size distribu- 
tion to an equilibrium power law. We chose a relatively 
steep value of a = —2.8 to decrease the computation 
time required for the size distribution evolution tests. 

As an initial test of the collisional evolution simu- 
lated by SMACK, we placed 1000 superparticles around 
a solar-mass star with semi-major axes uniformly dis- 
tributed between 90 AU and 110 AU. The superparti- 
cles were given eccentricities uniformly distributed be- 
tween 0 and e max = 0.2 and inclinations between 0 and 
i max = e-max/ 2. The other orbital elements were uni- 
formly distributed between 0 and 27 t. We distributed 
the planetesimals in each superparticle among 31 loga- 
rithmic size bins ranging from 1 mm in diameter to 1 m 
using a logarithmic increment of 0.1. In order to keep the 
average impact velocity constant to match the Dohnanyi 
(1969) scenario, we turned off the velocity corrections in 
SMACK by allowing SMACK to update the size distri- 
bution of each superparticle after an encounter without 
updating the superparticle trajectories. 

If the assumption of an infinite collisional cascade is 
relaxed and a cutoff is present in the small-sized end of 
the size distribution, a wave-like pattern emerges in the 
size distribution (Campo Bagatin et al. 1994). A real 
cutoff of this sort occurs in disks at the particle blowout 
size for the system, producing real wave patterns, but ar- 
tificial wave patterns can also appear as numerical noise 
in collisional models with artificial particle size cutoffs. 
We encountered this wave in our initial tests of the colli- 
sional algorithm in SMACK. Our simulated disk had an 
artificial cutoff in the size distribution at 1 mm. 

To remove this artificial wave, we use the method 
of Durda (1993). For each superparticle encounter, 
SMACK extrapolates the size distribution for each su- 
perparticle to 30 smaller size bins down to 1 /im. The 
number of collisions due to planetesimals in each of these 
smaller bins is calculated using equation (7) and is in- 
cluded in the particle loss and energy loss calculations 
(equations 15 and 10) for each superparticle. This small- 
particle extrapolation reduces the wave effect of the cut- 
off to negligible levels. 

We tested SMACK’s ability to reproduce the conver- 
gence to a Dohnanyi (1969) size distribution by running 
five simulations with different initial indices for the total 
size distributions for the planetesimals in the disk. For 
each simulation, set the vertical optical depth of the disk 
to 10~ 2 . We summed the size distribution over all super- 
particles and fit a power law to the total size distribution 
at each output tinrestep. Then we compared the evolu- 
tion of the power law indices over time for each of the 
five simulations. The results are shown in Fig. (1). As 
each collisional cascade evolved, the size distribution in- 
dex grew or shrank until it reached the Dohnanyi (1969) 
equilibrium index of —2.5. The size distribution index for 
each simulation converged to the Dohnanyi (1969) index 
within 2 x 10 6 yr. For comparison, the collision time for 
the smallest planetesimals is ~ 10 5 yr in this simulation. 

Using a size-dependent breaking strength (e.g. Gaspar 


et al. 2012b) or a size-dependent velocity distribution 
(see Pan & Schlichting 2012) can cause breaks in the 
steady-state distribution of a collisional cascade. We did 
not attempt to reproduce these phenomena in these ini- 
tial numerical tests. 



Fig. 1.- Incremental size distribution of the planetesimals in 
each of five different SMACK simulations. Each simulation had 
an initial power law size distribution of planetesimals between 1 
mm and 1 m, with power law indices ranging from po = —2.3 to 
po = —2.7. Each size distribution equillibrated to a power law with 
index p ps —2.5 within 2 X 10 6 yr (Dohnanyi 1969). 


3.2. Numerical Heating 

Simulations of bouncing collisions suffer from numer- 
ical heating associated with the finite size of the super- 
particles. The superparticles in SMACK do not sim- 
ply bounce, but the simulations of bouncing collisions 
can nonetheless inform us about the numerical noise in 
SMACK. For example, Lithwick & Chiang (2007) found 
that their grid-based simulations of bouncing collisions 
contained numerical heating on the scale of their grid 
size. Likewise, we expect that SMACK cannot model 
disks where the mean eccentricity of the planetesimals is 

< r s /r or the mean inclination of the planetesimals is 

< r s /r where r is the radius of the disk. 

We studied the eccentricity evolution of SMACK sim- 
ulations as a function of superparticle size to search for 
additional numerical heating effects, like any artificial 
viscous stirring associated with the finite superparticle 
size (see Goldreich & Tremaine 1978). We simulationed 
a planetless ring with radius 90-110 AU and optical depth 
5 x 10 -3 around a solar-mass star. We varied the super- 
particle radius between runs while varying the number of 
superparticles to keep the superparticle encounter time 
constant and measured how the mean eccentricity of the 
planetesimals evolved. We utilized the results of this test 
to choose the maximum superparticle size to use for the 
simulation described in Section 4. 

At first, we found that simulations using a large range 
of planetesimal sizes (1 mm - 1 m) cooled more slowly, 
and the eccentricity evolution curves never converged for 
any value of r s . This situation probably resulted from 
the coupling between planetesimals of various sizes. To 
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reduce this kind of numerical noise, we decided to restrict 
the size range of the planetesimals in each superparticle 
to 1 mm - 10 cm for all simulations that included velocity 
evolution. 

Fig. (2) shows the simulations with the reduced range 
of planetesimal sizes. The mean eccentricity of the su- 
perparticles in each simulation decreased rapidly at first 
due to inelastic collisions, then flattened and leveled off 
at a nonzero eccentricity. As we decreased the superpar- 
ticle radius, the eccentricity evolution curved converged 
to a single damping curve, in which the mean eccentricity 
dropped by a factor of 1/2 in 0.53 Myr, or approximately 
‘It co i , where t co i is the collision timescale for the smallest 
planetesimals. The damping curves converged at a su- 
perparticle radius of r s = 10 -13 AU, showing that this 
superparticle size probably suffices to mitigate numerical 
heating in our astrophysical simulations. 



Fig. 2. — The average eccentricity of the planetesimals in a ring 
with no planet, with radius 90-110 AU, over 10 Myr for a series 
of simulations that model the same disk using different superpar- 
ticle sizes. The damping rates converge for for superparticle radii 
S3 10 — 1,3 AU, suggesting that numerical heating is negligible for 
superparticles of this size. 


3.3. Numerical Viscosity 

The finite size of the superparticles also leads to dif- 
fusion of angular momentum, which can cause a narrow 
ring to spread on an artificially short timescale (e.g. Lith- 
wick & Chiang 2007). We tested our models for numer- 
ical viscous spreading by simulating the evolution of a 
very narrow (A r = 0.01 AU) planet-less ring at r = 10 
AU with superparticle radii of 10 -25 AU for 1 Myr. Fig. 
(3) shows the histogram of the number density of su- 
perparticles at several times during the simulation. The 
spreading time of this ring, the time it takes for the width 
of the ring to double, is T sprea d = 6.67 x 10 5 yr. 

Petit & Henon (1987) solved the diffusion equation for 
a collisional ring and found that the width of the ring A r 
at time t is related to the particle size (in our case, the 
superparticle size) r s by 



36kNr^t 

tperr 


1/3 


( 21 ) 



9.96 9.98 10.00 10.02 10.04 

Semi-Major Axis (AU) 


Fig. 3. — Histograms of superparticle semi-major axis at loga- 
rithmic time increments of a planetless ring at 10 AU with initial 
width 0.01 AU and superparticle radius 10 -2 ’ 5 AU. Numerical vis- 
cosity causes this extremely narrow ring to spread by a factor of 
two in 6.67 x 10 5 yr. 

where N is the number of particles, r is the radius of the 
ring, t per is the orbital period at that radius, and k is a 
dimensionless constant. This equation provides a good 
fit to the simulation shown in Fig. (3) when we choose 
k = 7.6 x 10 -4 . With this constant, we can use equation 
(21) to choose a superparticle size that yields acceptable 
levels of numerical viscosity. For example, if we want to 
limit the spreading of the ring in this simulation to 10% of 
its initial width in 10 7 yr, we must choose a superparticle 
size less than r s = 10 -412 AU. 


3.4. Conservation of Total Angular Momentum 

Collisions between planetesimals produce fragments 
smaller than the smallest size bin tracked by the super- 
particles. The mass of these dust particles is effectively 
lost to the superparticles, and carries away some angu- 
lar momentum. We measured the angular momentum of 
the superparticles and the mass lost to small fragments 
in the planet-less ring simulation described in Section 
3.2 to verify that these torques balanced and that total 
angular momentum was conserved by SMACK. 

We calculated the time derivative of the angular mo- 
mentum of the superparticles at each output timestep t 
as 


dL 

dt 


super particle 


y L^t + At) - y L^t) 

i i 

At 


( 22 ) 


where At is the size of an output timestep and Lift) is 
the magnitude of the angular momentum of superparticle 
i at time t, given by 


Lift) = nii(t)\vi(t)\\ri(t)\ sin 0j(t), (23) 


where r iff) and Vj(i) are the the position and veloc- 
ity vectors (in the heliocentric frame) of superparticle 
i, mift) is the total mass of superparticle i, and 9i(t) is 
the angle between rift) and Vj(i). 

We then set SMACK to output the mass lost in every 
superparticle encounter and the velocity of that mass. 
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We calculated the time derivative of the angular momen- 
tum of the small fragments lost as dust in the same way 
as we calculated the time derivative of the superparticle 
angular momentum: 


dC 


dt 


dust 


j k 


At 


(24) 


where Cj { t ) is the angular momentum of the small frag- 
ments produced in each of the j collisions during the 
output timestep t , and C k (t) the fragments produced in 
the k collisions during the previous output timestep. 

In Fig. 4 we plotted dL/dt of the superparticles, dC/dt 
of the lost mass, and the total for the two populations. 
The figure shows that most of the angular momentum 
lost by the superparticles is gained by the lost mass. The 
total change in angular momentum for the system varies 
stochastically over the simulation with a maximum vari- 
ation of 0.48% of the initial total angular momentum, Lq, 
per Myr. The systematic change in angular momentum 
is only 1.79 x 10~ 3 of L 0 over 10 Myr. We also calculated 
the change in the x, y, and z-components of the system’s 
angular momentum and found that the systemic change 
of each component was —3.13 x 10 _3 Lo, —4.46 x 10 _2 Lo, 
and 1.79 x 10~ 3 Lq over 10 Myr, respectively. 



Fig. 4. — The time-derivative of the angular momentum dL/dt 
of the disk over 10 Myr due to the superparticles and the mass 
lost to smaller fragments. The middle line shows the sum of dL/dt 
for the superparticles and lost mass. The total angular momentum 
decay rate fluctuates stochastically around zero, but the system 
lost only 0.179% of its initial angular momentum over the entire 
10 Myr simulation to numerical noise. 


3.5. Mass Loss Rates 


where Mq is the initial mass of the disk, t is a character- 
istic time at which the disk mass is half its initial mass, 
and C = 1/Mqt (Dominik & Decin 2003). 

We compared the evolution of the ring described in 
Section 3.2 to this analytic expression. Fig (5) shows the 
total mass of the disk during each simulation, with the 
full SMACK model and also for an identical model with 
velocity evolution turned off (i.e., the velocities of the 
superparticles were not updated). The mass loss curves 
for each simulation are in good agreement with the Do- 
minik & Decin (2003) prediction for t < t = 3.7 x 10 5 
yr, but the simulations lose mass at a faster rate than 
predicted for t > r. The mass loss curve for the full 
SMACK simulation begins to flatten at t ~ 10 6 yr as the 
mean eccentricity of the superparticles is damped (see 
Fig. (2)), decreasing the rate of fragmentation and slow- 
ing the mass loss. Lohne et al. (2008) also studied this 
problem with a model that included velocity evolution, 
but they did not report on this effect. More recently, 
Gaspar et al. (2012a) used a collisional model to demon- 
strate a similar mass loss damping to the damping seen 
in our simulations in Fig. (5). 



Fig. 5. — Total disk mass for a disk at a = 100 AU with width 
A a = 20 AU and initial maximum eccentricity e = 0.2 for a 
simulation without velocity evolution (dotted line) and with the 
full SMACK model (solid line). Both simulations are in good 
agreement with the Dominik & Decin (2003) analytic prediction 
of Mishit) = Mo/(l + £/ r ) (dashed line) for t < r = 1.4 x 10 5 yr, 
but the mass loss curves begin to diverge from from the analytic 
prediction for t > r. 

Rings of identical initial mass will have different char- 
acteristic times depending on their radii and the mean 
eccentricities of the planetesimals. Wyatt et al. (2007) 
analytically derived the following power law dependence 
of C on the radius of a ring: 


In the absence of velocity evolution, the rate of mass 
loss in a debris disk from planetesimal collisions is given 
by 

M disk (t) = -CM 2 disk , (25) 

and the total mass in a debris disk at time t is given by 

M dlsk (t) = (26) 


C oc r" 13/3 . (27) 

This power law is the product of three contributions: the 
r~ 1 / 2 dependence of the relative velocities in the ring, the 
r -5/6 dependence of the total number of projectiles above 
the minimum required disruptive size on the impact ve- 
locities, and an r~ 3 dependence of the density in the ring 
(Lohne et al. 2008). Lohne et al. (2008) studied mass loss 
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rates in debris disks and replicated this r“ 13 / 3 depen- 
dence in their simulations. In our simulations, the initial 
conditions of each system are determined by a user-input 
optical depth, so the density of the ring is independent 
of its size. This choice eliminates the r~ 3 factor for the 
ring density, so we predict a power law dependence on 
ring radius of: 

C oc r“ 4/3 . (28) 

We measured the mass loss in three simulated rings 
of radial width A r = 1 AU at five different radii from 
r = 1 to 3 AU, evolved for 10 4 yr. We turned off 
SMACK’s velocity evolution for these simulations by al- 
lowing SMACK to update the size distributions of the 
superparticles in each encounter without changing their 
velocities. We calculated C for each ring and fit a power 
law to the relationship between C and r for each ring 
(Fig. 6). The resulting dependence of C on radius in our 
simulated rings was C oc r -1 ' 24542 , in good agreement 
with the prediction of C oc r -4 / 3 . 



log Radius (AU) 

Fig. 6. — C = 1/Mot vs. disk radius for a disk with eccentricity 
e = 0.1, width A a = 1 AU, and no velocity evolution. The linear 
fit yields a power law with index -1.24542, close to the analytic 
prediction of —4/3. 

Wyatt et al. (2007) also derived the dependence of C 
on the mean eccentricity, e, of the ring planetesimals: 

C oc e 5 / 3 . (29) 

Lohne et al. (2008) also tested this mass loss rate depen- 
dence on eccentricity, but found a power law dependence 
of C oc e 9 / 4 . They argued that the analytical model must 
be incomplete. 

We measured the mass loss-eccentricity relationship in 
our models by simulating three rings at the same ra- 
dius (r = 1 AU) with different mean planetesimal eccen- 
tricities. Again, we turned off the velocity evolution in 
SMACK to compare the simulations with the predictions 
of Wyatt et al. (2007). After measuring the mass loss in 
each ring, we calculated C for each ring and fit a power 
law to the results (Fig. 7). The fit to the simulated data 
yields an index of 1.56767, in good agreement with the 


derived index of 5/3. However, at higher mean eccen- 
tricities, the radial excursions of the particles on their 
eccentric orbits can begin to dominate the width of the 
ring. At this point, increasing eccentricity will decrease 
the optical depth of the ring, which slows the mass loss 
rate dependence on eccentricity and flattens out the C 
vs. e relationship. Fig. (7) begins to show this effect in 
our simulations at higher eccentricities. 



log Eccentricity 

Fig. 7. — Disk eccentricity vs. C = 1/Mot for a disk with a = 1 
AU, width A a = 1 AU, and no velocity evolution. The linear 
fit yields a power law with index 1.56767, close to the analytic 
prediction of 5/3. 

4. A DISK CONTAINING A PLANET ON AN ECCENTRIC 

ORBIT 

To illustrate the application of SMACK to an inter- 
esting astrophysical system, we simulated the evolution 
of a debris ring containing a planet on an eccentric or- 
bit. We simulated a disk of planetesimals using 10,000 
superparticles with initial semi-major axes ranging from 
90-110 AU around a solar-mass star. We used a super- 
particle radius of r s = 10 -4 / 3 AU to maintain the same 
superparticle encounter time as the simulations in Fig. 
(2). Using equation (21) we can see that with this super- 
particle radius, numerical viscosity will cause the ring to 
spread by less than a factor of 10~ 4 in 10 Myr. 

The superparticle eccentricities and inclinations were 
uniformly distributed between 0 — 0.2 and 0 — 0.1, respec- 
tively. The remaining superparticle orbital elements U, 
w, and M, were distributed uniformly between 0 and 27 t. 
We inserted a planet with mass 8 Mj up at 25 AU with 
eccentricity 0.5 and zero inclination. The vertical opti- 
cal depth of the planetesimals at 100 AU from the star, 
Tdisk , w as set to 1 x 10 -2 . We evolved the system for 10 
Myr using SMACK. The simulation required only 4 wall 
clock hours on 48 cores on the NCCS Discover cluster. 

SMACK output the coordinates and size distributions 
of the superparticles every 10 4 yr. For each of four times 
during the simulation ( t = 0, 10 5 , 10 6 , and 10 7 yr), we 
used these coordinates and size distributions to synthe- 
size images of the disk. To synthesize each image, we 
combined together the superparticle data from 10 out- 
put timesteps to reduce the Poisson noise in the image 
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in the image, and projected the superparticles onto a 2D 
grid with resolution 2 AU to yield a face-on image. We 
then calculated each superparticlc’s contribution to the 
surface brightness in its pixel assuming each component 
planetesimal emits thermally as a spherical blackbody 
with 

luSP = fspTSpBviTsp), (30) 

where B v {Tsp ) is the value of the Planck function at 
Tgp, the temperature at the superparticle’s distance 
from the star. Then we summed over all the superparti- 
cles in each pixel to find the total surface brightness of 
that pixel, assuming a central star with solar luminosity. 
We simulated images at a wavelength of 850/im. The 
results of the are shown in Fig. (8). 

The first image in Fig. (8), in the upper left, shows the 
initial conditions of the ring when the planet is added. 
The ring is circular and centered on the star. In the next 
image at t = 10 5 yr, in the upper right, the planetesimal 
orbits are beginning to precess about the the planet’s 
forced eccentricity. At 10 5 yr the maximum vertical op- 
tical depth has dropped to r d i S k = 6.4 x 10~ 3 . Note that 
the ring is brightest at the top left corner of this frame, 
not at periapse to the right of the frame. This increased 
surface brightness in the top left corner probably arises 
from differential precession; planetesimals with smaller 
semi-major axes precess slightly faster than those in the 
outer part of the ring. 

By t = 10 6 yr (lower left frame), differential precession 
has sheared the disk out into a hint of a spiral structure, 
as described in Wyatt (1999). This spiral can be seen 
most clearly in Fig. (9), which uses a logarithmic bright- 
ness scale for a simulated image t = 5 x 10 5 yr. Then — 
in stark contrast to Wyatt (1999), which did not incor- 
porate collisions — collisions break up the spiral before 
it can wrap once more around the star, leaving an apse- 
aligned ring shown in the lower right, as described by 
Quillen (2006). The maximum vertical optical depth at 
t = 10 6 yr is Tdisk = 1-7 x 10~ 3 . The apse-aligned ring 
that remains at t = 10 7 yr (lower right frame) is now 
narrower, with a width of ~ 15 AU, because collisions re- 
move objects in orbits that are not approximately nested 
after precessing differentially for the lifetime of the sys- 
tem. At this final stage, the maximum vertical optical 
depth in the ring has dropped to Tdisk = 2.4 x 10 -4 . 

Inelastic collisions between particles can both excite 
and damp random velocities (Goldreich & Tremaine 
1978). Some authors have assumed that planetesimals in 
a debris disk have had their free eccentricities completely 
damped by collisions (Quillen 2006; Chiang et al. 2009) 
resulting in a ring of planetesimals on orbits with eccen- 
tricity equal to the forced eccentricity from the planet. 
Farther from the planet, the magnitude of the forced ec- 
centricity is given by 


\ e f 


orced \ 



(31) 


where e p is the planet’s eccentricity, by 2 (a) and b\, 2 (a) 
are Laplace coefficients, and a is the ratio of the planet’s 
semi-major axis a p to the ring’s semi-major axis a for 
a p < a (Murray & Dermott 1999). In the disk shown in 
Fig. (8), the forced eccentricity at a semi-major axis of 
100 AU is ~ 0.155. The mean eccentricity of the plan- 


etesimals in the t = 10 7 image in Fig. (8) is ~ 0.224, 
indicating that the free eccentricities have not been com- 
pletely damped, and the ring is not yet collisionally re- 
laxed. 

Besides secular perturbations and collisions, mean mo- 
tion resonances can also sculpt debris rings, clearing cen- 
tral holes (e.g. Mustill & Wyatt 2012) and creating a 
variety of azimuthal structures (e.g. Kuchner & Holman 
2003). However, the inner edge of the ring in our simula- 
tion is far outside the region of resonance overlap and an 
inspection of the time-averaged semi-major axes of the 
surviving superparticles near the end of the simulation 
(10 7 yr) suggests that only about 65 out of 1500 occupy 
the strongest mean motion resonances with the planet 
(7:1, 8:1, and 9:1) that intersect the ring. More sim- 
ulations will be needed to study the role of resonances 
in collisional rings, and the ability of SMACK to model 
such fine phase space structures using finite-sized super- 
particles. But for now, it appears that mean motion res- 
onances play no more than a minor role in the evolution 
of this system, and that secular and collisional evolution 
dominate the dynamics. 

To check if the erosion of the ring edges could be a 
numerical artifact, we ran the same simulation with four 
times as many superparticles. A histogram of the fi- 
nal semi-major axis distribution for the higher-fidelity 
simulation was nearly identical to final semi-major axis 
distribution of the simulation shown in Fig. (8), with 
differences consistent with Poisson noise. This compar- 
ison suggests that the narrowing of the ring illustrated 
in the lower right of Fig. (8) is real, caused by the com- 
bined effects of collisional damping and collisional (vis- 
cous) stirring of the differentially-precessing ring. 


5. CURRENT LIMITATIONS OF SMACK 

The main limitations of the current version of SMACK 
derive from using the superparticles to each represent 
many decades in planetesimal mass. Since the dominant 
size bins tend to be the smallest size bins and the “swap- 
ping” described in Section 2.2 tends to mix the size distri- 
butions among the superparticles, the current version of 
the code tends to make errors in the position and velocity 
distributions of larger bodies. For example, the current 
version of SMACK does not accurately model mass seg- 
regation and the velocity distributions of planetesimals 
of different sizes, as in the work of Pan & Schlichting 
( 2012 ). 

To investigate the level of mass segregation that 
SMACK can reproduce, we simulated the evolution of a 
planetless ring of superparticles with initial eccentricity 
of 0.1, as described in Section 3.2. We chose a super- 
particle radius of 0.1 AU to minimize numerical heating. 
We found that the eccentricity of the 1 mm planetesimals 
damped faster that the eccentricity of the 1 m planetes- 
imals by roughly a factor of 2, i.e., the planetesimals 
in SMACK do undergo some mass segregation, but not 
enough. 

In a real disk, the damping rate for planetesimals of 
size D by bodies of size s < D is 


1 N(s)s 3 
oc — 

7~damp 


(32) 


where Td am p is the damping timescale for the planetes- 


Collisions in Debris Disks 


11 



-150 -100 -50 0 50 100 150 



0.0 0.2 0.4 0.6 0.8 

Brightness (mjy/arcsec 2 ) 



-150 -100 -50 0 50 100 150 

Projected horizontal distance (AU) 



0.00 0.05 0.10 0.15 0.20 

Brightness (mJy/arcsec 2 ) 



-150 -100 -50 0 50 100 150 


I 


0.0 0.1 0.2 0.3 0.4 0.5 0.6 

Brightness (mjy/arcsec 2 ) 



-150 -100 -50 0 50 100 150 

Projected horizontal distance (AU) 



0.000 0.005 0.010 0.015 0.020 0.025 

Brightness (mJy/arcsec 2 ) 


Fig. 8. — Simulated 850 [im images of a disk with an 8 Mj up planet with semi-major axis 25 AU and eccentricity 0.5 at 0 yr, 10 5 yr, 10 6 
yr, and 10 7 yr. The white ellipse shows the planet’s orbit; the white x indicates the location of the star. The maximum optical depth of 
the ring at 10 7 yr (lower right frame) is 2.4 x 10“ 4 . 


imals and N(s) is the incremental size distribution of 
the disk with logarithmic size bins (Pan & Schlichting 
2012). Assuming a Dohnanyi (1969) size distribution 
of N(s) oc s -2 ’ 5 , the damping rate goes as s 0 5 D~ 1 . 
Collisions for which s « D will dominate, so the col- 
lision rate is proportional to D~ 05 . We therefore ex- 
pect the damping timescale for 1 m planetesimals to be 
larger than the timescale for 1 mm planetesimals by a 
factor of (10 3 ) 1 / 2 ~ 32, substantially larger than our 
SMACK model. In future versions of SMACK, we will 


explore using more superparticles and limiting the range 
of planetesimal sizes that each superparticle represents 
to improve SMACK’s ability to model mass segregation 
in disks. 

SMACK is based on the following assumptions: 

• No radiative forces: Besides causing issues with 
mass segregation, the use of superparticles also re- 
quires that all planetesimals in a superparticle be 
subject to the same forces. For example, radiation 
pressure and Poynting-Robertson drag forces both 
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indicates that SMACK did not update the velocities of the superparticles due to collisions, 
indicates that multiple simulations were run with different values of the given parameter. 
c Indicates that SMACK updated the superparticle velocities using Equations (5) and (6). 
d Indicates that SMACK updated the superparticles in some of the simulations but not others. 

TABLE 1 

Initial conditions for all simulations in Sections 3 and 4. 
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Fig. 9. — Simulated 850/im image of a disk with an 8Mj up planet 
with semi-major axis 25 AU and eccentricity 0.5 at 5 X 10 s yr, 
from the same simulation shown in Fig. (8). The white ellipse 
shows the planet’s orbit; the white x indicates the location of the 
star. The brightness scale in this image is logarithmic. The image 
shows a hint of spiral structure developing, before it is destroyed 
by collisions. 

vary with the size of the planetesimal. In this pa- 
per, we restrict our size distributions to planetesi- 
mals larger than 1 mm, for which radiative forces 
are not significant in a typical debris disk. 

• No interaction within a superparticle: SMACK also 
assumes that the planetesimals within a superpar- 
ticle do not interact. A cloud of planetesimals or- 
biting a star with the similar trajectories would 
not experience fragmenting collisions between plan- 
etesimals due to their low relative velocities, but 
they could still experience collisions with other out- 
comes, such as bouncing or sticking. Also, larger 
planetesimals in the cloud would gravitationally in- 
fluence the other bodies. SMACK does not sim- 
ulate gravity between planetesimals, so it cannot 
model grain growth or accretion by gravity. 

• No gravitational interaction between superparti- 
cles: Aside from preventing grain growth, the lack 
of gravitational interaction between large bodies in 
SMACK also prevents modeling of dynamical ef- 


fects such as viscous stirring. 

• Simplified crushing law: The crushing law used in 
this version of SMACK models only catastrophic 
collisions, with no sticking, cratering, or rebound- 
ing collisions. SMACK also assumes a single, mass- 
independent material strength, and calculates the 
size of the largest fragment independent of im- 
pact velocity. This simplified crushing law can 
reproduce a size distribution evolution to a basic 
Dohnanyi (1969) equiliubrium. Future versions of 
SMACK will include more complex crushing laws 
based on experimental results and SPH modeling 
to investigate the effects of crushing law parame- 
ters on the size distribution evolution of a disk. 

6. SUMMARY AND FUTURE WORK 

SMACK, a new collisional module for the REBOUND 
N-body integrator, uses a superparticle method to sim- 
ulate the evolution of the dynamics and size distribu- 
tion of a debris disk experiencing fragmenting collisions. 
SMACK models the interaction between the disk and em- 
bedded planets in 3-D for planetesimals with a range of 
sizes and can model various disk asymmetries including 
eccentric rings and warps. When run on parallel CPUs, 
SMACK can simulate the evolution of disks older than 
10 7 yr in a feasible number of hours. 

We showed that SMACK, with its simple “swapping” 
algorithm, can reproduce a variety of fundamental re- 
sults in debris disk physics. With the velocity evolu- 
tion turned off (i.e., by not updating superparticle veloc- 
ities) , we used SMACK to reproduce the evolution of the 
size distribution of planetesimals to a power-law with 
index —2.5, using incremental logarithmic size bins, as 
Dohnanyi (1969) derived analytically. We verified that 
SMACK conserves total angular momentum during dust- 
producing collisions and showed that we could mitigate 
the effects of numerical heating and numerical viscosity 
by constraining the size of the superparticles. We repro- 
duced the mass loss function for a debris disk derived by 
Dominik & Decin (2003) and demonstrated that the rate 
of mass loss varies with disk radius and planetesimals 
mean eccentricity as roughly r -4 / 3 and e 5 / 3 , as Wyatt 
(2008) derived. 

After performing this battery of tests, we used SMACK 
to simulate the evolution of a disk containing a central 
planet on an eccentric orbit (Figure 8). The simula- 
tion shows how planetesimal collisions, together with the 
planet’s gravity, converted a circular ring of planetesi- 
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mals (width 20 AU), centered on the star, into a nar- 
rower (6 AU) ring, offset from the star. It illustrates the 
combined effects of secular precession, collisional damp- 
ing and collisional (viscous) stirring. Previous models of 
a planet on an eccentric orbit interacting with a debris 
disk have assumed either no particular collisional evolu- 
tion of the planetesimal orbits (Wyatt 1999; Wilner et al. 
2002; Quillen & Thorndike 2002; Moran et al. 2004), or 
complete collisional damping (Quillen 2006; Chiang et al. 
2009). Our model hints at the wide range of intermediate 
possibilities that are likely important in observed debris 
disks; in our simulation, the planetesimals were still not 
completely collisionally damped after 10' yr, roughly the 
age of f3 Pictoris. 

Boley et al. (2012) have suggested that a pair of planets 
might be needed to create the narrow planetesimal ring 
around Fomalhaut. But no second planet was required 
in our simulation. Of course, the scenario we simulated 
was ultimately an artificial one: a planet on an eccentric 
orbit suddenly introduced into a disk of planetesimals in 
a circular ring. More simulations, with a wider range 
of initial conditions and planet-formation scenarios are 
needed to unravel the physics of the Fomalhaut system. 

The primary limitation of the SMACK models pre- 
sented in this paper at that they underestimate the mass 
segregation rate, as described in section (5). We plan to 
try improving the way SMACK models planetesimal-size- 
dependent effects by limiting the range of planetesimal 
masses that each superparticle represents. We hope to 
extend the large end of the size distributions to 1 km and 
beyond. 


We anticipate several other future improvements. The 
current version of SMACK uses a size-independent 
breaking strength for the planetesimals, which we will re- 
place with a size-dependent strength to explore its effects 
on the steady-state size distribution of a collisional cas- 
cade. We will also expand the set of collisional outcomes 
to include cratering and bouncing collisions. Aggrega- 
tion and planetesimals growth could also conceivably be 
simulated with SMACK. 

We also intend to implement a dynamic and adaptive 
domain decompositioning for REBQUND to enable us to 
run simulations on even more CPUs. This will allow us 
to improve the running time of REBOUND and SMACK 
and to model even larger systems. 

With the flexibility of the REBOUND platform, 
SMACK can model a wide range of debris disk scenar- 
ios: disks with planets on eccentric orbits (e.g. Section 
4), disks with planets on inclined orbits like the /3 Pic- 
toris system (Lagrange et al. 2009, 2010), even multiple 
planet systems and migrating planets. Millimeter and 
sub- millimeter observations with ALMA are beginning to 
probe populations of larger dust grains and parent bod- 
ies in debris disks at high angular resolution (e.g. Boley 
et al. 2012; MacGregor et al. 2013). We hope to make 
SMACK a workhorse for interpreting ALMA images of 
debris disks. 

We thank the NASA High-end Computing Program for 
granting us time on the Discover cluster. This research 
was supported in part by NASA Planetary Geology and 
Geophysics Program, grant no. 11-PGG11-0032. 
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